Efficient segmentation of piecewise smooth images

ABSTRACT

A fast and robust segmentation model for piecewise smooth images is provided. Local statistics in an energy formulation are provided as a functional. The shape gradient of this new functional gives a contour evolution controlled by local averaging of image intensities inside and outside the contour. Fast computation is realized by expressing terms as the result of convolutions implemented via recursive filters. Results are similar to the general Mumford-Shah model but realized faster without having to solve a Poisson partial differential equation at each iteration. Examples are provided. A system to implement segmentation methods is also provided.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional Application No.60/857,295 filed Nov. 7, 2006, which is incorporated herein byreference.

BACKGROUND OF THE INVENTION

The present invention relates to image segmentation. In particular itrelates to efficient segmentation of images with smooth intensityregions.

The extraction of piecewise smooth regions from an image is of greatinterest in different domains, and still remains a challenging task. Forexample, this is very useful in medical imaging where organs orstructures of interest are often characterized by smooth intensityregions. This problem has been formulated as the minimization of anenergy by Mumford and Shah in: D. Mumford and J. Shah. Optimalapproximations by piecewise smooth functions and associated variationalproblems. Comm. Pure Appl. Math., 42:577-685, 1989 as:

E ^(MS)(u,Γ)=μ²∫_(Ω)(u ₀ −u)² dx+∫ _(Ω\Γ) |∇u| ² dx+ν|Γ|  (1)

where u is the piecewise smooth function, Γ is the interface betweensmooth regions, and u₀ is the original image. The interpretation of thethree terms is straightforward: the first one is the usual mean-squaredata term; the second one means that one wants to extract smoothregions; the third one means that one want to extract regions withsmooth boundaries.

The minimizer of this so-called Mumford-Shah functional gives a boundarythat separates the image domain in smooth regions. A very interestingproperty of this approach is that it solves two common image-processingtasks simultaneously: image denoising and image segmentation. However,finding the minimizer is not straightforward and remains an issue. Beingnon-convex, the functional is most of the time minimized using gradientdescent techniques which are subject to local minima. For example, in“L. Vese and T. Chan. A multiphase level set framework for imagesegmentation using the mumford and shah model. International Journal ofComputer Vision, 50:271-293, 2002” and “A. Tsai, A. Jr. Yezzi, and A. S.Willsky. Curve evolution implementation of the mumford-shah functionalfor image segmentation, denoising, interpolation, and magnification.IEEE Transactions on Image Processing, 10(8):1169-1186, August 2001”,the optimization process alternates between the evolution of one or twolevel set functions such as provided in “S. Osher and J. Sethian. Frontspropagation with curvature dependent speed: Algorithms based onHamilton-Jacobi formulations. J. of Comp. Phys., 79:12-49, 1988”, “J.Sethian. Level Set Methods and Fast Marching Methods: EvolvingInterfaces in Computational Geometry, Fluid Mechanics, Computer Vision,and Materials Sciences. Cambridge Monograph on Applied and ComputationalMathematics. Cambridge University Press, 1999” and “A. Dervieux and F.Thomasset. A finite element method for the simulation of Raleigh-Taylorinstability. Springer Lect. Notes in Math., 771:145-158, 1979” and theresolution of Poison partial differential equations. This process iscomputational extensive and requires a very good initialization to avoidbeing stuck into local minima.

To relax this problem, one can consider a restriction of E^(MS) topiecewise constant functions. Let Ω_(i) be the open subsets delimited byΓ, the piecewise constant Mumford-Shah functional writes:

$\begin{matrix}{{E_{0}^{MS}(\Gamma)} = {{\sum\limits_{i}{\int_{\Omega_{i}}{\left( {u_{0} - {{mean}_{\Omega_{i}}\left( u_{0} \right)}} \right)^{2}{x}}}} + {v{\Gamma }}}} & (2)\end{matrix}$

This functional was shown in “D. Mumford and J. Shah. Optimalapproximations by piecewise smooth functions and associated variationalproblems. Comm. Pure Appl. Math., 42:577-685, 1989” to be a limitfunctional of expression (1) as 0. A level set implementation of thisfunctional known as the Chan and Vese model was proposed in “L. Vese andT. Chan. A multiphase level set framework for image segmentation usingthe mumford and shah model. International Journal of Computer Vision,50:271-293, 2002”. While this simplified functional is easier tominimize, it also makes a very strong assumption on the image byassuming implicitly a Gaussian intensity distribution for each regionΩ_(i) as described in “M. Rousson and R. Deriche. A variationalframework for active and adaptative segmentation of vector valuedimages. In Proc. IEEE Workshop on Motion and Video Computing, pages56-62, Orlando, Fla., December 2002”. Other papers model thisdistribution with Gaussian mixtures such as in “N. Paragios and R.Deriche. Geodesic active regions and level set methods for supervisedtexture segmentation. The International Journal of Computer Vision,46(3):223-247, 2002” and “0. Juan, R. Keriven, and G. Postelnicu.Stochastic motion and the level set method in computer vision:Stochastic active contours. The International Journal of ComputerVision, 69(1):7-25, 2006”, or with nonparametric distributions such as“J. Kim, J. Fisher, A. Yezzi, M. Cetin, and A. Willsky. Nonparametricmethods for image segmentation using information theory and curveevolution, in IEEE International Conference on Image Processing, pages797-800, September 2002”, but they all make the assumption of a globaldistribution over each region. In many real images, these globalintensity models are not valid. This is often the case in medicalimaging, especially in MR images where an intensity bias can beobserved.

Several approaches are available to overcome the limitation of globaltechniques. One is to consider image gradients by fitting the contour toimage discontinuities. This is generally referred to as edge-basedmethods, and it is the basis of the Geodesic Active Contours asdescribed in “V. Caselles, R. Kimmel, and G. Sapiro. Geodesic activecontours. International Journal of Computer Vision, Vol 22:61-79, 1997”and “S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, and A. Yezzi.Gradient flows and geometric active contour models. In Proceedings ofthe 5^(th) International Conference on Computer Vision, pages 810-815,Boston, Mass., June 1995. IEEE Computer Society Press”. A thirdfunctional was also introduced in the seminal work of Mumford and Shah:“D. Mumford and J. Shah. Optimal approximations by piecewise smoothfunctions and associated variational problems. Comm. Pure Appl. Math.,42: 577-685, 1989”. This functional is the integral along I′ of ageneralized Finsler metric and leads indeed to the first geodesic activecontour (before “V. Caselles, R. Kimmel, and G. Sapiro. Geodesic activecontours. International Journal of Computer Vision, Vol 22:61-79, 1997”and “S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, and A. Yezzi.Gradient flows and geometric active contour models. In Proceedings ofthe 5^(th) International Conference on Computer Vision, pages 810-815,Boston, Mass., June 1995. IEEE Computer Society Press”.).

Edge-based methods are also well-known for their high sensitivity tonoise and for the presence of local minima in the optimization asdescribed in D. Cremers, M. Rousson, and R. Deriche. A review ofstatistical approaches to level set segmentation: integrating color,texture, motion and shape. International Journal of Computer Vision,2006. Another alternative was briefly discussed in “L. Vese. Multiphaseobject detection and image segmentation, in S. Osher and N. Paragios,editors, Geometric Level Set Methods in Imaging, Vision and Graphics,pages 175-194. Springer Verlag, 2003” where the function u₀ of theMumford-Shah functional is restricted to a linear function of thespatial location x: u₀(x)=a.x+b. Even though this last one is promising,it is still restricted to very particular spatial distributions of theintensity.

Accordingly methods and systems for extracting piecewise smooth regionsbut with improved performance are required.

SUMMARY OF THE INVENTION

In accordance with one aspect of the present invention a novel methodand system are provided for an improved segmentation of an image objectfrom a background.

In accordance with another aspect of the present invention a method isprovided for segmentation of an image u₀ comprising the steps ofcreating a plurality of smooth regions Ω_(i) delimited by a boundary Γ;approximating the image by a piecewise smooth function u_(σ) inside eachregion Ω_(i); and evaluating a functional.

In accordance with a further aspect of the present invention a method isprovided wherein the piecewise smooth function is expressed as:

u σ  ( x , Ω i ) = ∫ Ω i  g σ  ( x - y )  u 0  ( y )   y ∫  g σ ( x - y )   y ,  wherein   g σ  ( v ) = 1 2  πσ  exp  ( - v 22  σ 2 )

is a Gaussian kernel.

In accordance with another aspect of the present invention a method isprovided wherein the functional is expressed asE(F)=μ²∫_(S)(u₀−u_(σ)(Γ))²dx+ν|Γ|.

In accordance with a further aspect of the present invention a method isprovided wherein the functional is expressed as

E  ( Γ ) = μ 2  ∑ i  ∫  ( u 0 - u σ  ( ) ) 2   x + v   Γ  .

In accordance with another aspect of the present invention a method isprovided wherein u_(σ) represents an overall piecewise smoothapproximation expressed as

${u_{\sigma}\left( {x,\Gamma} \right)} = {\sum\limits_{i}{\chi_{i}{{u_{\sigma}\left( {x,\Omega_{i}} \right)}.}}}$

In accordance with a further aspect of the present invention a method isprovided wherein a performance of segmentation can be tuned by selectinga value of σ.

In accordance with another aspect of the present invention a method isprovided wherein the segmentation is a bi-partitioning, separating aregion Ω from a complementary region Ω.

In accordance with a further aspect of the present invention a method isprovided further comprising: evolving the boundary Γ between region Ωand region Ω according to an expression:

${\frac{\partial\Gamma}{\partial t}(x)} = {\left\lbrack {\mu^{2}\left( {\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\overset{\_}{\Omega}} \right)}} \right)^{2} - \left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)^{2} - {q_{\sigma}\left( {x,\overset{\_}{\Omega}} \right)} + q_{\sigma}} \right)} \right) + {v\; \kappa \text{]}{{N(x)}.}}}$

In accordance with another aspect of the present invention a method isprovided further comprising: implementing the evolving of the boundary Γwith a level set representation.

In accordance with a further aspect of the present invention a method isprovided further comprising computing steps of the evolving of theboundary Γ with a recursive filter.

In accordance with another aspect of the present invention a system isprovided which can perform and execute the steps of the methods hereprovided as aspects of the present invention.

DESCRIPTION OF THE DRAWINGS

FIG. 1 shows an original image and smoothed versions of that image.

FIG. 2 shows an original image with segmentation initialization and asegmented image.

FIG. 3 shows an image with a segmentation initialization.

FIG. 4. shows segmentation of the image of FIG. 3 in accordance with anaspect of the present invention.

FIG. 5. shows another segmentation of the image of FIG. 3 in accordancewith an aspect of the present invention.

FIG. 6 shows yet another segmentation of the image of FIG. 3 inaccordance with an aspect of the present invention.

FIG. 7 shows an image with a segmentation initialization.

FIG. 8 shows a segmentation of the image of FIG. 7.

FIG. 9 shows a segmentation of the image of FIG. 7 in accordance with anaspect of the present invention.

FIG. 10 illustrates a segmentation of an object in a medical image.

FIG. 11 illustrates a segmentation of an object in a medical image inaccordance with an aspect of the present invention.

FIG. 12 illustrates a computer system that is used to perform the stepsdescribed herein in accordance with another aspect of the presentinvention.

DESCRIPTION OF A PREFERRED EMBODIMENT

As an aspect of the present invention a general approach for extractingpiecewise smooth regions is disclosed. Instead of minimizing thedistance between the intensity and the average intensity of the regionlike in expression (2), the distance between the intensity and localaveraging inside the region is minimized. This gives a model able toapproximate piecewise smooth functions like the original Mumford-Shahfunctional provided in expression (1), but with a complexity closer tothat of the piecewise constant model.

The model here provided as another aspect of the present invention willbe explained in detail and how it can be linked to the Mumford-Shahmodel. The minimized energy as well as its derivative using the shapegradient are expressed. The level set method is used to compute theevolution of the interface, and each term of the derivative is expressedas the result of a convolution as a further aspect of the presentinvention. The importance of the fast recursive filter is brieflyexplained. Some results on synthetic and real data will be provided andcompared with the piecewise smooth and piecewise constant Mumford-ShahModel.

Piecewise Smooth Approximation Through Gaussian Convolutions

The Mumford-Shah model approximates the image by a piecewise smoothfunction by penalizing the high gradient of an “ideal” cartoon image (uin expression (1)). Making an estimation of this image at the same timeas the segmentation makes the functional difficult to minimize, andcomputationally expensive, as it is the solution of a poisson equation.

Herein the problem is approached differently by fixing the cartoon imageto a smoothing of the image intensity inside each subset Ω_(i) for the 2Dimensional case, the smooth function in Ω_(i) is then defined as,

u σ  ( x , Ω i ) = ∫ Ω i  g σ  ( x - y )  u 0  ( y )   y ∫  g σ ( x - y )   y ,  wherein   g σ  ( v ) = 1 2  πσ  exp  ( - v 22  σ 2 ) ( 3 )

The denominator of this expression is a normalization factor which isimportant for voxels that are close to the border of Ω_(i), i.e. whenthe Gaussian kernel does not overlap completely with Ω_(i). This can bealso interpreted as a local weighted averaging of the image intensityaround the voxel x inside the region Ω_(i). This is demonstrated inFIG. 1. It shows the importance of the denominator in the formulation ofsmooth regions. Image 101 in FIG. 1 is the original image; image 102 andimage 103 in FIG. 1 are smoothed approximations of the image 101 wherein102 and 103 are the piecewise smooth images estimated by the model whichis an aspect of the present invention, for a given position of theboundary. Estimate 102 is an approximation without the normalizationfactor and 103 is an approximation with the normalization factor.

Let χ_(i) be a characteristic function of Ω_(i) such that χ_(i)(x)=1 ifx∈

and 0 otherwise. One can express the overall piecewise smoothapproximation of the image as

$\begin{matrix}{{u_{\sigma}\left( {x,\Gamma} \right)} = {\sum\limits_{i}{\chi_{i}{u_{\sigma}\left( {x,\Omega_{i}} \right)}}}} & (4)\end{matrix}$

With this approximation, u is a piecewise smooth function that is givenanalytically with respect to the boundary Γ, and one does not needanymore the regularization term on u present in expression (1). Thisleads to a new functional:

E  ( Γ ) = μ 2  ∫ S  ( u 0 - u σ  ( Γ ) ) 2   x + v   Γ  = μ 2 ∑ i  ∫  ( u 0 - u σ  ( ) ) 2   x + v   Γ  ( 5 )

Interestingly, when the variance σ of the Gaussian goes to infinity,this functional becomes equivalent to the piecewise constant model. Thislimit model has become very popular in its level set formulation(Chan-Vese model as described in L. Vese and T. Chan. A multiphase levelset framework for image segmentation using the mumford and shah model.International Journal of Computer Vision, 50:271-293, 2002) because itperforms very well for regions that are characterized by quite differentglobal means. However it is not able to discriminate regions with nearlythe same global intensity distributions as is shown in FIG. 2. FIG. 2 isan example of an image that does not suit the Chan-Vese model. Image 201of FIG. 2 shows the initialization and image 202 the convergence. With adifferent choice of σ, the model here presented becomes more local andcan segment a wider set of images where image regions differ only intheir local intensity distributions. Hence, tuning the parameter σpermits to control the locality of intensity statistics inside eachregion.

Although the model here presented has no restriction on the number ofregions to segment, in the following the focus is on the bi-partitioningcase to make the explanations simpler. In particular, this allows torepresent the boundary Γ with a single level set function. Severalextensions using multiple level set functions have been proposed tosegment an arbitrary number of regions as described in “L. Vese and T.Chan. A multiphase level set framework for image segmentation using themumford and shah model. International Journal of Computer Vision,50:271-293, 2002”.

Energy Minimization

In the case of bi-partitioning, the contour Γ separates a region Ω fromits complementary Ω, and the energy (5) becomes:

$\begin{matrix}{{E(\Gamma)} = {{\mu^{2}{\sum\limits_{D = {\{{\Omega,\overset{\_}{\Omega}}\}}}{\int_{D}{\left( {{u_{0}(x)} - \frac{\int_{D}{{g_{\sigma}\left( {x - y} \right)}{u_{0}(y)}{y}}}{\int_{D}{{g_{\sigma}\left( {x - y} \right)}{y}}}} \right){x}}}}} + {v{\Gamma }}}} & (6)\end{matrix}$

To minimize this energy, one can use the shape gradient tools developedin “G. Aubert, M. Barlaud, 0. Faugeras, and S. Jehan-Besson. Imagesegmentation using active contours: Calculus of variations or shapegradients? SIAM Journal of Applied Mathematics, 63(6):2128-2154, 2003”.The detailed derivation is presented in a later section of thedisclosure. It leads to the following evolution on the boundary Γ:

$\begin{matrix}{{{{\frac{\partial\Gamma}{\partial t}(x)} = {\left\lbrack {{\mu^{2}\left( {\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\overset{\_}{\Omega}} \right)}} \right)^{2} - \left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)^{2} - {q_{\sigma}\left( {x,\overset{\_}{\Omega}} \right)} + {q_{\sigma}\left( {x,\Omega} \right)}} \right)} + {v\; \kappa}} \right\rbrack {N(x)}}},\mspace{20mu} {with}}{{q_{\sigma}\left( {x,\Omega} \right)} = {\int_{\Omega}{\frac{2\left( {{u_{0}(y)} - {u_{\sigma}\left( {u,\Omega} \right)}} \right)\left( {{u_{0}(x)} - {u_{\sigma}\left( {y,\Omega} \right)}} \right){g_{\sigma}\left( {y - x} \right)}}{\int_{\Omega}{{g_{\sigma}\left( {y - z} \right)}{z}}}{{y}.}}}}} & (7)\end{matrix}$

Where N(x) is the outward normal vector to F at the point x. The firsttwo terms of this evolution equation are similar to the ones that can befound in the usual piecewise smooth and constant Mumford-Shah cases.Their interpretation is quite straightforward: the contour will locallymove to include the current image voxel in the region it is the mostsimilar to. The other terms are unique to the formulation provided hereas an aspect of the present invention and come from the analyticalexpression of the piecewise smooth image as a function of the boundary.

Level Set Implementation

Any curve representation can be used to implement the evolutiondescribed in expression (7). Here it is presented how to do this with alevel set representation. In particular, this allows giving animplementation that is valid in any dimension. Any curve representationcan be used to implement the evolution described in expression (7). Hereit is presented how to do it with a level set representation. Inparticular, this allows to give an implementation that is valid in anydimension.

Let φ be the signed distance function to the boundary Γ, positive in Ωand negative in Ω. H_(α) is here introduced as the regularized versionsof the Heaviside functions. Equation (7) becomes:

$\begin{matrix}\left. {\frac{\partial\varphi}{\partial t} = {{{{\nabla\varphi}}\left\lbrack {{\mu^{2}\left( {\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\overset{\_}{\Omega}} \right)}} \right)^{2} - \left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)^{2} - {q_{\sigma}\left( {x,\overset{\_}{\Omega}} \right)}} \right)} + {q_{\sigma}\left( {x,\Omega} \right)}} \right)} + {v\; \kappa}}} \right\rbrack & (8)\end{matrix}$

where all four different terms u_(σ)(φ),ū_(σ)(φ),q_(σ)(φ) and q _(σ)(φ)can be computed with convolutions by the Gaussian kernel g_(σ):

$\begin{matrix}\left\{ \begin{matrix}{{u_{\sigma}(\varphi)} = \frac{g_{\sigma}*{H_{\alpha}(\varphi)}u_{0}}{g_{\sigma}*{H_{\alpha}(\varphi)}}} \\{{{\overset{\_}{u}}_{\sigma}(\varphi)} = \frac{g_{\sigma}*\left( {1 - {H_{\alpha}(\varphi)}} \right)u_{0}}{g_{\sigma}*\left( {1 - {H_{\alpha}(\varphi)}} \right)}} \\{{q_{\sigma}(\varphi)} = {{u_{0}\left( {g_{\sigma}*\frac{2\left( {u_{0} - u_{\sigma}} \right)}{g_{\sigma}*{H_{\alpha}(\varphi)}}} \right)} - {g_{\sigma}*\frac{2\left( {u_{0} - u_{\sigma}} \right)u_{\sigma}}{g_{\sigma}*{H_{\alpha}(\varphi)}}}}} \\{{{\overset{\_}{q}}_{\sigma}(\varphi)} = {{u_{0}\left( {g_{\sigma}*\frac{2\left( {u_{0} - {\overset{\_}{u}}_{\sigma}} \right)}{g_{\sigma}*\left( {1 - {H_{\alpha}(\varphi)}} \right)}} \right)} - {g_{\sigma}*\frac{2\left( {u_{0} - {\overset{\_}{u}}_{\sigma}} \right){\overset{\_}{u}}_{\sigma}}{g_{\sigma}*\left( {1 - {H_{\alpha}(\varphi)}} \right)}}}}\end{matrix} \right. & (9)\end{matrix}$

Each one of these terms need to be updated at each evolution of thelevel set. Even though these expressions seem complicated, theirestimation is quite straightforward since they are the results ofconvolutions by a Gaussian kernel (more details are provided in a latersection). This is a good advantage because it can be implemented veryeffectively with a recursive filter as described in “R. Deriche. UsingCanny's criteria to derive a recursively implemented optimal edgedetector. The International Journal of Computer Vision, 1(2):167-187,May 1987”. Hence, for a d dimensional image with N voxels, thecomplexity of each convolution is in O(d. N), only six convolutions areneeded to compute the four terms. If one compares this complexity to thepiecewise constant case, where the means inside each region also need tobe recomputed at each iteration, it is of the same order.

Results and Comparisons

The model provided as an aspect of the present invention is demonstratedby applying it on several images. FIG. 3 shows the initialization of animage. FIGS. 4-6 show the role of the variance σ in evolution: a smallvariance extracts thinner details, whereas a large variance suits wellfor large homogeneous shape. As pointed out earlier, when σ goes toinfinity, the functional become equivalent to the Chan-Vese model.However, for small variances, the initial contour has to be close toedges to evolve, but are able to extract much thinner details. Actually,the model behaves like the geodesic active contour model without anyballoon force: it drives the front toward edges in the image, and makesit evolve via a mean curvature flow in homogeneous regions. Moregenerally, the evolution of the front “follows” the edges that cross theinitial contour. The model is thus very dependent on the initialization.

FIGS. 7-9 demonstrate an aspect of the present invention. FIG. 7 showsthe initialization of the segmentation of a leaf object. The object inFIG. 7 has two distinct regions characterized by the same globalstatistics. FIG. 8 shows the segmentation using the Chan-Vese model. InFIG. 8, one can point out the limitations of the Chan-Vese model thatseparates the white regions from the black, also in the background, anddoes not extract accurately the leaf in the image. FIG. 9 shows asegmentation of the leaf object in accordance with an aspect of thepresent invention. One can see in FIG. 9 that the segmentation modelhere provided model behaves quite well. One can clearly see that thefront “follows” the edges that crossed it in its initial shape.

In FIGS. 10 and 11, both the model which is an aspect of the presentinvention in FIG. 11 and the Chan-Vese model in FIG. 10 are applied inorder to extract a 2D liver from biased anatomical MRI. As the Chan-Vesemodel is not robust to bias, the liver is not correctly extracted, andthe front “leak” in the part of the image where intensity is close tothe global mean inside the front as can be seen in FIG. 10. One can seein FIG. 11 that the model which is an aspect of the present inventionsbehaves quite well, as the liver is correctly delimited at the locallevel. One just has to make the initial front cross the edges of theliver, and make the front evolve, thus following the edges of the liver,and finally extract it almost completely.

In conclusion, a new model was presented for extracting smooth regionsfrom image data. This model is based on the Mumford-Shah functional, butis formulated in a simpler and more efficient way. A new functional wasprovided, and shown that it was linked to the Chan-Vese model, byrepresenting regions as a local average instead of a global mean. Oneinteresting point is that the minimization of this functional can becomputed in a very fast way, thanks to the Deriche recursive filter.Promising results were provided on 2D synthetic and real data.

Derivation of Energy Equation

In an earlier section, an energy was defined that is composed of twodomain integrals. Here, a detailed derivation is provided using shapegradients. The following energy is considered:

$\begin{matrix}{{{E\left( {\Omega,\overset{\_}{\Omega}} \right)} = {{\int_{\Omega}{{f\left( {x,\Omega} \right)}{x}}} + {\int_{\overset{\_}{\Omega}}{{f\left( {x,\overset{\_}{\Omega}} \right)}{x}}}}},{where}} & \; \\{{{f\left( {x,\Omega} \right)} = {\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)^{2} = \left( {{u_{0}(x)} - \frac{G_{1}\left( {x,\Omega} \right)}{G\; 2\left( {x,\Omega} \right)}} \right)^{2}}},{and}} & (10) \\\left\{ \begin{matrix}{{G_{1}\left( {x,\Omega} \right)} = {{\int_{\Omega}{{H_{1}\left( {x,y} \right)}{y}\mspace{14mu} {and}\mspace{14mu} {H_{1}\left( {x,y} \right)}}} = {{u_{0}(y)}{g_{\sigma}\left( {x - y} \right)}}}} \\{{G_{2}\left( {x,\Omega} \right)} = {{\int_{\Omega}{{H_{2}\left( {x,y} \right)}{y}\mspace{14mu} {and}\mspace{14mu} {H_{2}\left( {x,y} \right)}}} = {g_{\sigma}\left( {x - y} \right)}}}\end{matrix} \right. & (11)\end{matrix}$

Recalling the main theorem of shape gradient method presented in “G.Aubert, M. Barlaud, 0. Faugeras, and S. Jehan-Besson. Image segmentationusing active contours: Calculus of variations or shape gradients? SIAMJournal of Applied Mathematics, 63(6):2128-2154, 2003”:

Theorem 1. The Gâteaux derivative of a functional J(Ω)=g(x,Ω)dx in thedirection of a vectorfield V is:

J′(Ω), V

=∫ _(Ω) g _(s)(x,Ω,V)dx−∫ _(∂Ω) g(x,Ω)(V(x).N(x))da(x),

where g_(x)(x,Ω,V) is the shape derivative of g(x,Ω) in the direction ofV; ∂Ω is the boundary of Ω; N is the unit inward normal to ∂Ω; and daits area element.

In the present case, for the bi-partitioning problem, one has,

${\langle{{E^{\prime}\left( {\Omega,\overset{\_}{\Omega}} \right)},V}\rangle} = {{\int_{\Omega}{{f_{s}\left( {x,\Omega,C} \right)}{x}}} + {\int_{\overset{\_}{\Omega}}{{f_{s}\left( {x,\overset{\_}{\Omega},V} \right)}{x}}} - {\int_{\Gamma}{\left( {{f\left( {x,\Omega} \right)} - {f\left( {x,\overset{\_}{\Omega}} \right)}} \right)\left( {{V(x)} \cdot {N(x)}} \right){{a(x)}}}}}$  withf_(s)(x, Ω, V) = f_(G₁)(x, Ω, G₁, G₂)⟨G₁^(′)(x, Ω) ⋅ V⟩ + f_(G₂)(x, Ω, G₁, G₂)⟨G₂^(′)(x, Ω) ⋅ V⟩,

where f_(G) ₁ and f_(G2) denote the partial derivative of equation (10)with respect to G₁ and G₂. They can be expressed as:

$\begin{matrix}\left\{ \begin{matrix}{{f_{G_{1}}\left( {x,\Omega,G_{1},G_{2}} \right)} = {{- \frac{2}{G_{2}\left( {x,\Omega} \right)}}\left( {{I(x)} - \frac{G_{1}\left( {x,\Omega} \right)}{G_{2}\left( {x,\Omega} \right)}} \right)}} \\{{f_{G_{2}}\left( {x,\Omega,G_{1},G_{2}} \right)} = {{- \frac{2{G_{1}\left( {x,\Omega} \right)}}{{G_{2}\left( {x,\Omega} \right)}^{2}}}\left( {{I(x)} - \frac{G_{1}\left( {x,\Omega} \right)}{G_{2}\left( {x,\Omega} \right)}} \right)}}\end{matrix} \right. & (12)\end{matrix}$

and by using Theorem 1,

$\begin{matrix}\left\{ \begin{matrix}{{\langle{{G_{1}^{\prime}\left( {x,\Omega} \right)} \cdot V}\rangle} = {{\int_{\Omega}{{H_{1s}\left( {x,y,V} \right)}{y}}} - {\int_{\Gamma}{{H_{1}\left( {x,y} \right)}\left( {{V(y)} \cdot {N(y)}} \right){{a(y)}}}}}} \\{{\langle{{G_{2}^{\prime}\left( {x,\Omega} \right)} \cdot V}\rangle} = {{\int_{\Omega}{{H_{2s}\left( {x,y,V} \right)}{y}}} - {\int_{\Gamma}{{H_{2}\left( {x,y} \right)}\left( {{V(y)} \cdot {N(y)}} \right){{a(y)}}}}}}\end{matrix} \right. & (13)\end{matrix}$

Since H₁ and H₂ do not depend on Ω, one obtains H_(1s)=0 and H_(2s)=0.Putting all the terms together one finds:

$\begin{matrix}\begin{matrix}{{f_{s}\left( {x,\Omega,V} \right)} = {{f_{G_{1}}{\langle{G_{1}^{\prime} \cdot V}\rangle}} + {f_{G_{2}}{\langle{G_{2}^{\prime} \cdot V}\rangle}}}} \\{= \frac{2\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)}{\int_{\Omega}{{g_{\sigma}\left( {x - y} \right)}{y}}}} \\{{\int_{\Gamma}{{g_{\sigma}\left( {x - y} \right)}\begin{pmatrix}{{u_{0}(y)} -} \\{u_{\sigma}\left( {x,\Omega} \right)}\end{pmatrix}\begin{pmatrix}{{V(y)} \cdot} \\{N(y)}\end{pmatrix}{{a(y)}}}}}\end{matrix} & \begin{matrix}(14) \\(15)\end{matrix}\end{matrix}$

and at last one obtains by changing the order of integration,

$\begin{matrix}{\begin{matrix}{{\langle{{E^{\prime}\left( {\Omega,\overset{\_}{\Omega}} \right)},V}\rangle} = {{\int_{\Omega}{{f_{s}\left( {x,\Omega,V} \right)}{x}}} + {\int_{\overset{\_}{\Omega}}{{f_{s}\left( {x,\overset{\_}{\Omega},V} \right)}{x}}} -}} \\{{\int_{\Gamma}\left( {{f\left( {y,\Omega} \right)} - {\left( {f\left( {y,\overset{\_}{\Omega}} \right)} \right)\left( {{V(y)} \cdot {N(y)}} \right){{a(y)}}}} \right.}} \\{= {\int_{\Gamma}\left( {{q\left( {y,\Omega} \right)} - {q\left( {y,\overset{\_}{\Omega}} \right)} - \begin{pmatrix}{{u_{0}(y)} -} \\{u_{\sigma}\left( {y,\Omega} \right)}\end{pmatrix}^{2} +} \right.}} \\{\left. \left( {{u_{0}(y)} - {u_{\sigma}\left( {y,\overset{\_}{\Omega}} \right)}} \right)^{2} \right)\left( {{V(y)} \cdot {N(y)}} \right){{a(y)}}}\end{matrix}{with}{{q\left( {y,\Omega} \right)} = {\int_{\Omega}{\frac{2\begin{pmatrix}{{u_{0}(x)} -} \\{u_{\sigma}\left( {x,\Omega} \right)}\end{pmatrix}{g_{\sigma}\left( {x - y} \right)}\begin{pmatrix}{{u_{0}(y)} -} \\{u_{\sigma}\left( {x,\Omega} \right)}\end{pmatrix}}{\int_{\Omega}{{g_{\sigma}\left( {x - z} \right)}{z}}}{x}}}}} & (16)\end{matrix}$

One finally gets the following gradient descent:

$\begin{matrix}{\left. {{\frac{\partial\Gamma}{\partial t} (x)} = \left\lbrack {\left( {{q\left( {x,\Omega} \right)} - {q\left( {x,\overset{\_}{\Omega}} \right)} - \left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)^{2} + \mspace{445mu} \left( {{{u_{0}( x)} - u_{\sigma}}, \overset{\_}{\Omega}} \right)} \right)}^{2} \right)} \right\rbrack  {N( x)}} & (17)\end{matrix}$

Implementation

Each integral term in the gradient descent can be seen as a convolutionby a kernel of variance σ. Therefore, the computation of the evolutioncan be done in a very fast way in two separate steps:

first: make several convolutions via a fast recursive filter (refderiche)

next: compute the speed of each point in the narrow-band by using thepreviously computed blurred images.

One has,

$\begin{matrix}{{{u\left( {x,\Omega} \right)} = {\frac{\int_{\Omega}{{g_{\sigma}\left( {x - y} \right)}{u_{0}(y)}{y}}}{\int_{\Omega}{{g_{\sigma}\left( {x - y} \right)}{y}}} = \frac{\left( {g_{\sigma}*u_{0}} \right)_{\Omega}(x)}{\left( {g_{\sigma}*1} \right)_{\Omega}(x)}}}{and}\begin{matrix}{{q\left( {y,\Omega} \right)} = {\int_{\Omega}\frac{2\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right){g_{\sigma}\left( {x - y} \right)}\left( {{u_{0}(y)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)}{\int_{\Omega}{{g_{\sigma}\left( {x - z} \right)}{z}}}}} \\{= {{{u_{0}(y)}{\int_{\Omega}{\frac{2\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)}{\left( {g_{\sigma}*1} \right)_{\Omega}(x)}{g_{\sigma}\left( {x - y} \right)}{x}}}} -}} \\{{\int_{\Omega}{\frac{2\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right){u_{\sigma}\left( {x,\Omega} \right)}}{\left( {g_{\sigma}*1} \right)_{\Omega}(x)}{g_{\sigma}\left( {x - y} \right)}{x}}}} \\{= {{u_{0}(y)}\left( {{g_{\sigma}*q_{1}}_{\Omega}{{(y) - \left( {g_{\sigma}*q_{2}} \right)}_{\Omega}(y)}} \right.}}\end{matrix}{with}} & (18) \\{{{q\; 1\left( {x,\Omega} \right)} = \frac{2\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)}{\left( {g_{\sigma}*1} \right)_{\Omega}(x)}}{and}{{q\; 2\left( {x,\Omega} \right)} = \frac{2\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right){u_{\sigma}\left( {x,\Omega} \right)}}{\left( {g_{\sigma}*1} \right)_{\Omega}(x)}}} & (19)\end{matrix}$

The domain convolution can be computed by using the Heaviside H_(σ) ofthe level-set function:

$\begin{matrix}\left\{ \begin{matrix}{{\left( {g_{\sigma}*f} \right)_{\Omega}(x)} = {\left( {g_{\sigma}*H_{\sigma}f} \right)(x)}} \\{{\left( {g_{\sigma}*f} \right)_{\overset{\_}{\Omega}}(x)} = {\left( {g_{\sigma}*\left( {1 - H_{\sigma}} \right)f} \right)(x)}}\end{matrix} \right. & (20)\end{matrix}$

System

The segmentation methods that are aspects of the present invention canbe executed by a system as shown in FIG. 12. The system is provided withdata 1201 representing image data. An instruction set or program 1202executing the methods of the present invention is provided and combinedwith the data in a processor 1203, which can process the instructions of1202 applied to the data 1201. A segmented image can be outputted andfor instance displayed on a display 1204. The processor can be dedicatedhardware. However, the processor can also be a CPU or any othercomputing device that can execute the instructions of 1202. An inputdevice 1205 like a mouse, or track-ball or other input device allows auser to select an initial object to be segmented and to start thesegmenting process. Accordingly the system as shown in FIG. 12 providesa system for efficient image segmentation using methods disclosedherein.

The following references are generally descriptive of the background ofthe present invention and are hereby incorporated herein by reference:[1]. G. Aubert, M. Barlaud, 0. Faugeras, and S. Jehan-Besson. Imagesegmentation using active contours: Calculus of variations or shapegradients? SIAM Journal of Applied Mathematics, 63(6):2128-2154, 2003.[2]. V. Caselles, R. Kimmel, and G. Sapiro. Geodesic active contours.International Journal of Computer Vision, Vol 22:61-79, 1997. [3]. D.Cremers, M. Rousson, and R. Deriche. A review of statistical approachesto level set segmentation: integrating color, texture, motion and shape.International Journal of Computer Vision, 2006. [4]. R. Deriche. UsingCanny's criteria to derive a recursively implemented optimal edgedetector. The International Journal of Computer Vision, 1(2):167-187,May 1987. [5]. A. Dervieux and F. Thomasset. A finite element method forthe simulation of Raleigh-Taylor instability. Springer Lect. Notes inMath., 771:145-158, 1979. [6]. 0. Juan, R. Keriven, and G. Postelnicu.Stochastic motion and the level set method in computer vision:Stochastic active contours. The International Journal of ComputerVision, 69(1):7-25, 2006. [7]. S. Kichenassamy, A. Kumar, P. Olver, A.Tannenbaum, and A. Yezzi. Gradient flows and geometric active contourmodels. In Proceedings of the 5^(th) International Conference onComputer Vision, pages 810-815, Boston, Mass., June 1995. IEEE ComputerSociety Press. [8]. J. Kim, J. Fisher, A. Yezzi, M. Cetin, and A.Willsky. Nonparametric methods for image segmentation using informationtheory and curve evolution. In IEEE International Conference on ImageProcessing, pages 797-800, September 2002. [9]. D. Mumford and J. Shah.Optimal approximations by piecewise smooth functions and associatedvariational problems. Comm. Pure Appl. Math., 42: 577-685, 1989. [10].S. Osher and J. Sethian. Fronts propagation with curvature dependentspeed: Algorithms based on Hamilton-Jacobi formulations. J. of Comp.Phys., 79:12-49, 1988. [11]. N. Paragios and R. Deriche. Geodesic activeregions and level set methods for supervised texture segmentation. TheInternational Journal of Computer Vision, 46(3):223-247, 2002. [12]. M.Rousson and R. Deriche. A variational framework for active andadaptative segmentation of vector valued images. In Proc. IEEE Workshopon Motion and Video Computing, pages 56-62, Orlando, Fla., December2002. [13]. J. Sethian. Level Set Methods and Fast Marching Methods:Evolving Interfaces in Computational Geometry, Fluid Mechanics, ComputerVision, and Materials Sciences. Cambridge Monograph on Applied andComputational Mathematics. Cambridge University Press, 1999. [14]. A.Tsai, A. Jr. Yezzi, and A. S. Willsky. Curve evolution implementation ofthe mumford-shah functional for image segmentation, denoising,interpolation, and magnification. IEEE Transactions on Image Processing,10(8):1169-1186, August 2001. [15]. L. Vese. Multiphase object detectionand image segmentation. In S. Osher and N. Paragios, editors, GeometricLevel Set Methods in Imaging, Vision and Graphics, pages 175-194.Springer Verlag, 2003. [16]. L. Vese and T. Chan. A multiphase level setframework for image segmentation using the mumford and shah model.International Journal of Computer Vision, 50:271-293, 2002.

The terms voxel and voxels used herein are also intended to mean pixelor pixels in the appropriate situations.

While there have been shown, described and pointed out fundamental novelfeatures of the invention as applied to preferred embodiments thereof,it will be understood that various omissions and substitutions andchanges in the form and details of the device illustrated and in itsoperation may be made by those skilled in the art without departing fromthe spirit of the invention. It is the intention, therefore, to belimited only as indicated by the scope of the claims appended hereto.

1. A method of segmentation of an image u₀ comprising: creating aplurality of smooth regions Ω_(i) delimited by a boundary Γ;approximating the image by a piecewise smooth function u_(σ) inside eachregion Ω_(i); and evaluating a functional.
 2. The method as claimed inclaim 1, wherein the functional is evaluated in a level setimplementation.
 3. The method as claimed in claim 1, wherein thepiecewise smooth function is expressed as: u σ  ( x , Ω i ) = ∫ Ω i  gσ  ( x - y )  u 0  ( y )   y ∫  g σ  ( x - y )   y ,  wherein  g σ  ( y ) = 1 2  πσ  exp  ( - v 2 2  σ 2 ) is a Gaussiankernel.
 4. The method as claimed in claim 1, wherein the functional isexpressed asE(Γ)=μ²˜_(S)(u₀ −u _(σ)(σ))² dx+ν|Γ|.
 5. The method as claimed in claim1, wherein the functional is expressed as E  ( Γ ) = μ 2  ∑ i  ∫  (u 0 - u σ  ( ) ) 2   x + v   Γ  .
 6. The method as claimed inclaim 1, wherein u_(σ) represents an overall piecewise smoothapproximation expressed as${u_{\sigma}\left( {x,\Gamma} \right)} = {\sum\limits_{i}{\chi_{i}{{u_{\sigma}\left( {x,\Omega_{i}} \right)}.}}}$7. The method as claimed in claim 3, wherein a performance ofsegmentation can be tuned by selecting a value of σ.
 8. The method asclaimed in claim 1, wherein the segmentation is a bi-partitioning,separating a region Ω from a complementary region Ω.
 9. The method asclaimed in claim 8, further comprising: evolving the boundary Γ betweenregion Ω and region Ω according to an expression${\frac{\partial\Gamma}{\partial t}(x)} = {\left\lbrack {{\mu^{2}\left( {\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\overset{\_}{\Omega}} \right)}} \right)^{2} - \left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)^{2} - {q_{\sigma}\left( {x,\overset{\_}{\Omega}} \right)} + {q_{\sigma}\left( {x,\Omega} \right)}} \right)} + {v\; \kappa}} \right\rbrack {{N(x)}.}}$10. The method as claimed in claim 9, further comprising: implementingthe evolving of the boundary Γ with a level set representation.
 11. Themethod as claimed in claim 10, further comprising: computing steps ofthe evolving of the boundary Γ with a recursive filter.
 12. A system ofsegmentation of an image u₀ comprising: a processor; software operableon the processor to: creating a plurality of smooth regions Ω_(i)delimited by a boundary Γ; approximating the image by a piecewise smoothfunction u_(σ) inside each region Ω_(i); and evaluating a functional.13. The system as claimed in claim 12, wherein the functional isevaluated in a level set implementation.
 14. The system as claimed inclaim 12, wherein the piecewise smooth function is expressed as: u σ  (x , Ω i ) = ∫ Ω i  g σ  ( x - y )  u 0  ( y )   y ∫  g σ  ( x -y )   y ,  wherein   g σ  ( v ) = 1 2  πσ  exp  ( - v 2 2  σ 2) is a Gaussian kernel.
 15. The system as claimed in claim 12, whereinthe functional is expressed asE(Γ)=μ²∫_(S)(u ₀ −u _(σ)(γ))² dx+ν|Γ|.
 16. The system as claimed inclaim 12, wherein the functional is expressed as E  ( Γ ) = μ 2  ∑ i ∫  ( u 0 - u 0  ( ) ) 2   x + v   Γ  .
 17. The system as claimedin claim 12, wherein u_(σ) represents an overall piecewise smoothapproximation expressed as${u_{\sigma}\left( {x,\Gamma} \right)} = {\sum\limits_{i}{\chi_{i}{{u_{\sigma}\left( {x,\Omega_{i}} \right)}.}}}$18. The system as claimed in claim 14, wherein a performance ofsegmentation can be tuned by selecting a value of σ.
 19. The system asclaimed in claim 12, wherein the segmentation is a bi-partitioning,separating a region Ω from a complementary region Ω.
 20. The system asclaimed in claim 19, further comprising: evolving the boundary Γ betweenregion Ω and region Ω according to an expression${\frac{\partial\Gamma}{\partial t}(x)} = {\left\lbrack {{\mu^{2}\left( {\left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\overset{\_}{\Omega}} \right)}} \right)^{2} - \left( {{u_{0}(x)} - {u_{\sigma}\left( {x,\Omega} \right)}} \right)^{2} - {q_{\sigma}\left( {x,\overset{\_}{\Omega}} \right)} + {q_{\sigma}\left( {x,\Omega} \right)}} \right)} + {v\; \kappa}} \right\rbrack {{N(x)}.}}$21. The system as claimed in claim 20, further comprising: implementingthe evolving of the boundary Γ with a level set representation.
 22. Thesystem as claimed in claim 21, further comprising: computing steps ofthe evolving of the boundary Γ with a recursive filter.